########################################
#Replication code for "Regions at Risk"#
#by Sebastian Schutte (2015)           #
########################################

#How to run this code:
#1. Create a new folder at some suitable location. Place the replication data into this folder. 

#2. Change the following line to your replication folder with the replication data. The code produces multiple figures, 
#   but the tables are only produced as latex code in the R output. You can save them with write.csv().

setwd('/home/basti/Desktop/psrm_rep2/new_repdata/psrm_repdata')

#3. Execute the following block of 'library' commands

library (maptools)
library (spatstat)
library (raster)
library (plotKML)
library (xtable)

#4. Execute the following block of helper functions to do the simulations and cross validations

###############
#Aux functions#
###############
to_ppp <- function(acled,cshape){
  acled  <- as(acled,"SpatialPoints")
  acled  <- as(acled,"ppp")
  acled  <- unique(acled)
  cshape <- as(cshape,"SpatialPolygons")
  cshape <- as(cshape,"owin")
  ppp  <- ppp(acled$x, acled$y,window=cshape)
  return(ppp)
}
to_ppm <- function(country,polygon,n_predictors){
  im.pop <- as.im(as(crop(pop, polygon),'SpatialGridDataFrame'))
  im.lsg <- as.im(as(crop(lsg, polygon),'SpatialGridDataFrame'))
  im.acc <- as.im(as(crop(acc, polygon),'SpatialGridDataFrame'))
  im.veg <- as.im(as(crop(veg, polygon),'SpatialGridDataFrame'))
  im.nrd <- as.im(as(crop(nrd, polygon),'SpatialGridDataFrame'))
  im.brd <- as.im(as(crop(brd, polygon),'SpatialGridDataFrame'))
  
  if(n_predictors == 1){
    return (ppm (country, ~  pop ,
                 covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg)))
  } 
  if(n_predictors == 2){
    return (ppm (country, ~  pop + lsg,
                 covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg)))
  }
  if(n_predictors == 3){
    return (ppm (country, ~  pop + lsg + acc,
                 covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg)))
  }  
  if(n_predictors == 4){
    return (ppm (country, ~  pop + lsg + acc + brd,
                 covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg)))
  }
  if(n_predictors == 5){
    return (ppm (country, ~  pop + lsg + acc + brd + nrd,
                 covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg)))
  }
  if(n_predictors == 6){
    return (ppm (country, ~  pop + lsg + acc + brd + nrd + veg,
                 covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg)))
  }
}
to_ppm_star <- function(country,polygon){
  im.pop <- as.im(as(crop(pop, polygon),'SpatialGridDataFrame'))
  im.lsg <- as.im(as(crop(lsg, polygon),'SpatialGridDataFrame'))
  im.acc <- as.im(as(crop(acc, polygon),'SpatialGridDataFrame'))
  im.veg <- as.im(as(crop(veg, polygon),'SpatialGridDataFrame'))
  im.nrd <- as.im(as(crop(nrd, polygon),'SpatialGridDataFrame'))
  im.brd <- as.im(as(crop(brd, polygon),'SpatialGridDataFrame'))
  return (ppm (country, ~  pop + nrd  + brd + acc,
               covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg)))
}
norm_dens <- function(current_im){
  max_value <- max(current_im$v,na.rm=T)
  if(max_value > 0){
    current_im$v[is.na(current_im$v)] <- 0
    current_im$v <- current_im$v / max_value
    current_im$v[current_im$v == 0] <- NA
  }
  return(current_im)
}
#################
#Simulate Events#
#################
simulate_events <- function(country,MLE,country_name,n_sim,csr){
  sim_runs <- c()
  for (i in 1:n_sim){
    if (csr == TRUE){
      sim      <- density(rpoint(country$n,win=country))
    } else {
      sim      <- density(rmh(model=MLE,start=list(n.start=country$n), control=list(p=1,nrep=1000,nverb=5),new.coef=MLE$coef))    
    }
    sim_runs <- append(sim_runs,as.matrix(sim))
    gc()
  }
  result <- c(rep(0,length(sim_runs)/n_sim))
  offset <- length(result)
  agg    <- 0
  for (i in 1:length(result)){
    for (j in 0:(n_sim-1)){
      agg <- agg + sim_runs[i+(j*offset)]
    } 
    result[i] <- agg / n_sim
    agg <- 0
  }
  result_plot <- matrix(nrow=sqrt(length(result)),ncol=sqrt(length(result)))
  for (i in 1:length(result_plot[1,])){
    for (j in 0:length(result_plot[,1]-1)){
      result_plot[i,j] <- result[i + (length(result_plot[,1]) * j)]
    }
  } 
  cntry <- density(country)
  cntry <- as.matrix(cntry)
  
  kmllist <- list()
  kmlnames <- list()
  kmllist[[1]] <- raster(density(country))
  kmlnames[[1]] <- paste(country_name,"-emp",sep="")
  kmlPlot(kmllist, kmlnames, paste(country_name,"-emp.kml",sep=""))

  kmllist[[1]] <- raster(im(result_plot))
  kmlnames[[1]] <- paste(country_name,"-sim",sep="")
  kmlPlot(kmllist, kmlnames, paste(country_name,"-sim.kml",sep=""))
}

####################
#Make general model#
####################

make_general_model <- function(country_list,country_fit_list){
  k <- length(country_list)
  general_fit <- country_fit_list[[1]]
  coefs <- country_fit_list[[1]]$coef
  for (i in 1:k){
    coefs <- cbind(coefs,country_fit_list[[i]]$coef)  
  }
  for (i in 1:length(general_fit$coef)){  
    general_fit$coef[i] <- mean(coefs[i,])   
  }
  return (general_fit$coef)
}

############################
#Get differences in density#
############################

get_density_difference <- function(country,country_fit,country_name,n_sim,csr){
  sim_runs <- c()
  for (i in 1:n_sim){
    if (csr == TRUE){
      sim      <- density(rpoint(country$n,win=country)) 
    } else { 
      sim      <- density(rmh(model=country_fit,start=list(n.start=country$n), control=list(p=1,nrep=1000,nverb=5),new.coef=country_fit$coef))
    }    
    sim_runs <- append(sim_runs,as.matrix(sim))
    gc()
    cat(i,"/",n_sim,country_name,"\n")
  }
  sim_runs <- as.vector(sim_runs)
  sim_results     <- c(rep(0,length(sim_runs)/n_sim))
  sim_sd_results  <- c(rep(0,length(sim_runs)/n_sim))
  offset <- length(sim_results)
  
  for (i in 1:length(sim_results)){
    agg <- c()
    for (j in 0:(n_sim-1)){
      agg <- c(agg,sim_runs[i+(j*offset)])
    }
    agg[is.na(agg)] <- 0
    sim_results[i]     <- mean(agg)
    sim_sd_results[i]  <- sd(agg)
  }
  mae <- c(rep(0,offset))
  emp_result <- as.vector(as.matrix(density(country)))
  emp_result[is.na(emp_result)] <- 0
  sim_results <- sim_results / max(sim_results)   
  emp_result <- emp_result / max(emp_result)
  mae <- mean(abs(emp_result - sim_results)) 
  return(mae)
}

##################
#Cross validation#
##################
k_fold_cross_plot <- function (countries_list,country_fit_list,country_names,csr,n_sim){
  k <- length(countries_list)
  result_matrix <- as.data.frame(matrix(nrow=k,ncol=2))#country, csr, model
  names(result_matrix) <- c("country","cv")
  for(i in 1:k){
    cat("#############################",country_names[[i]],"\n")
    #Split the sample in training and testing part
    predict_country         <- countries_list[[i]]
    predict_fit             <- country_fit_list[[i]]
    training_countries       <- countries_list[-i]
    training_fits            <- country_fit_list[-i]
    traning_coefs <- make_general_model(training_countries,training_fits)
    predict_fit$coef[1:length(predict_fit$coef)]   <- traning_coefs[1:length(predict_fit$coef)]
    if (csr == TRUE){
      simulate_events(predict_country,predict_fit,paste(country_names[[i]],"_CSR"),n_sim,csr)
    } else {
      simulate_events(predict_country,predict_fit,paste(country_names[[i]],"_CV"),n_sim,csr)
    }
    result_matrix$country[i] <- country_names[[i]]
    result_matrix$cv[i] <- get_density_difference(predict_country,predict_fit,country_names[[i]],n_sim,csr)          
    gc()
  }
  return (result_matrix)
}

is_predict_plot <- function(countries_list,country_fit_list,country_names,csr,n_sim){
  k <- length(countries_list)
  result_matrix <- as.data.frame(matrix(nrow=k,ncol=2))#country, csr, model
  names(result_matrix) <- c("country","cv")
  for(i in 1:k){
    result_matrix$country[i] <- country_names[[i]]
    result_matrix$cv[i] <- get_density_difference(countries_list[[i]],country_fit_list[[i]],conflict_name_list[[i]],n_sim,csr)
    gc()
  }
  return(result_matrix)
}

kmlPlot <- function(datalist, namelist, kmlfilename){
  k <- 1
  plotKML.env(silent = TRUE, colour_scale_numeric = "", kmz = FALSE)
  kml_open(file.name = kmlfilename,overwrite = TRUE)
  while(k <= length(datalist)){
    proj4string(datalist[[k]]) <- CRS("+proj=longlat +datum=WGS84")
    kml_layer.Raster(datalist[[k]],plot.legend=FALSE,colour_scale = SAGA_pal[[1]], alpha = 1,raster_name = paste(namelist[[k]],"-country",".png",sep=""))
    k <- k + 1
  }
  kml_close(kmlfilename)
}

get_var_gain <- function(country_fit, conflict_data, is, n_sim){
  im.pop <- as.im(as(crop(pop,(as((conflict_data$window),'SpatialPolygons'))),'SpatialGridDataFrame'))
  im.lsg <- as.im(as(crop(lsg,(as((conflict_data$window),'SpatialPolygons'))),'SpatialGridDataFrame'))
  im.veg <- as.im(as(crop(veg,(as((conflict_data$window),'SpatialPolygons'))),'SpatialGridDataFrame'))
  im.acc <- as.im(as(crop(acc,(as((conflict_data$window),'SpatialPolygons'))),'SpatialGridDataFrame'))
  im.nrd <- as.im(as(crop(nrd,(as((conflict_data$window),'SpatialPolygons'))),'SpatialGridDataFrame'))
  im.brd <- as.im(as(crop(brd,(as((conflict_data$window),'SpatialPolygons'))),'SpatialGridDataFrame'))
  
  permutations <- get_binary_permutations6()
  permutations <- cbind(permutations,NA)
  colnames(permutations) <- c("pop","lsg","acc","brd","nrd","veg","cv")
  for (i in 1:length(permutations[,1])){
    cat(i,"/",length(permutations[,1]),"\n")
    callstring <- ""
    for (j in 1:6){
      if (permutations[i,j] == 1){
        callstring <- paste(callstring,colnames(permutations)[j])
      }
    }
    callstring <- paste(strsplit(trim(callstring),split=" ")[[1]],collapse =" + ")
    callstring <- paste(" ~",callstring)#,",covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg)")
    #specify trend
    m <- ppm(conflict_data, trend=as.formula(callstring), 
             covariates = list (pop = im.pop, lsg = im.lsg, acc = im.acc, nrd = im.nrd, brd = im.brd,veg = im.veg))
    if (is == TRUE){
      permutations[i,7] <- is_predict_plot(list(conflict_data),list(m),list("none"),FALSE,n_sim)$cv
    } else {
      permutations[i,7] <- k_fold_cross_plot(list(conflict_data),list(m),list("none"),FALSE,n_sim)$cv
    }
  }
  return(data.frame(permutations))
}

#Helper function to generate binary permutations. Borrowed from http://stackoverflow.com/questions/5671149/permute-all-unique-enumerations-of-a-vector-in-r
get_binary_permutations6 <- function(){
  p <- c(1,1,1,1,1,1)
  p <- rbind(p,
             uniqueperm2(c(1,1,1,1,1,0)),
             uniqueperm2(c(1,1,1,1,0,0)),
             uniqueperm2(c(1,1,1,0,0,0)),
             uniqueperm2(c(1,1,0,0,0,0)),
             uniqueperm2(c(1,0,0,0,0,0)))            
  return(p)
}
uniqueperm2 <- function(d) {
  dat <- factor(d)
  N <- length(dat)
  n <- tabulate(dat)
  ng <- length(n)
  if(ng==1) return(d)
  a <- N-c(0,cumsum(n))[-(ng+1)]
  foo <- lapply(1:ng, function(i) matrix(combn(a[i],n[i]),nrow=n[i]))
  out <- matrix(NA, nrow=N, ncol=prod(sapply(foo, ncol)))
  xxx <- c(0,cumsum(sapply(foo, nrow)))
  xxx <- cbind(xxx[-length(xxx)]+1, xxx[-1])
  miss <- matrix(1:N,ncol=1)
  for(i in seq_len(length(foo)-1)) {
    l1 <- foo[[i]]
    nn <- ncol(miss)
    miss <- matrix(rep(miss, ncol(l1)), nrow=nrow(miss))
    k <- (rep(0:(ncol(miss)-1), each=nrow(l1)))*nrow(miss) + 
      l1[,rep(1:ncol(l1), each=nn)]
    out[xxx[i,1]:xxx[i,2],] <- matrix(miss[k], ncol=ncol(miss))
    miss <- matrix(miss[-k], ncol=ncol(miss))
  }
  k <- length(foo)
  out[xxx[k,1]:xxx[k,2],] <- miss
  out <- out[rank(as.numeric(dat), ties="first"),]
  foo <- cbind(as.vector(out), as.vector(col(out)))
  out[foo] <- d
  t(out)
}

#5. Set the variable below to a suitable value. To replicate the results from the manuscript, leave it at 100. This will 
#   result in very long run times, however, because 100 point patterns will be simulated for each model and each country 
#   to assess predictive gain. To just generate some quick results that will differ slightly from the published ones
#   set the value to something like 3 or 5. 

n_simulations <- 100

#6. Now run the below code that uses the above helper functions to generate figures and tables 
#   from the manuscript. 

#############
#Read cshape#
#############
cshapes <- readShapeSpatial("cshapes.shp", proj4string=CRS("+proj=longlat +ellps=WGS84"))
date <- as.Date("2002-1-1")
startdate <- as.Date(paste(cshapes$COWSYEAR, cshapes$COWSMONTH, cshapes$COWSDAY, sep="-"))
enddate <- as.Date(paste(cshapes$COWEYEAR, cshapes$COWEMONTH, cshapes$COWEDAY, sep="-"))
cshp2000 <- cshapes[startdate < date & enddate >= date,]
africa <- readShapePoly("africa.shp",proj4string=CRS("+proj=longlat +ellps=WGS84"))
prediction_countries <- subset(cshp2000, is.element(GWCODE,africa$gwnum) | is.element(GWCODE,c(645,652,770,700,630,645,652,660,666,663,670,698,678)))

#################
#Read covariates#
#################
pop    <- raster ("population.tif",proj4string = CRS("+proj=longlat +ellps=WGS84"))
acc    <- raster ("access.tif",proj4string = CRS("+proj=longlat +ellps=WGS84"))
lsg    <- raster ("capdist.tif",proj4string = CRS("+proj=longlat +ellps=WGS84"))
nrd    <- raster ("wealth.tif",proj4string = CRS("+proj=longlat +ellps=WGS84"))
brd    <- raster ("borderdist.tif",proj4string = CRS("+proj=longlat +ellps=WGS84"))
veg    <- raster ("landcover.tif",proj4string = CRS("+proj=longlat +ellps=WGS84"))

gc()

##########
#Read ged#
##########
gedshp <- readShapePoints("ucdp-ged-15-export.shp",proj4string=CRS("+proj=longlat +ellps=WGS84"))
ins_data <- read.csv("idata.csv")
ged <- subset(gedshp, is.element(RELID, ins_data$relid))
relevant_countries <- subset(cshp2000, is.element(GWCODE,unique(ins_data$gwno)))
#Remove Djibouti (see manuscript)
relevant_countries <- subset(relevant_countries, ISONAME != 'Djibouti')

#I set the random seed for replication purposes:
set.seed(12345)

######################
#Generate model lists#
######################

conflict_data_list  <- list()
conflict_model_list_1 <- list()
conflict_model_list_2 <- list()
conflict_model_list_3 <- list()
conflict_model_list_4 <- list()
conflict_model_list_5 <- list()
conflict_model_list_6 <- list()
conflict_model_list_star <- list()

conflict_name_list  <- list()

k <- 1
while(k <= length(relevant_countries$GWCODE)){
	c_id <- relevant_countries$GWCODE[k]
  g_tmp <- ged[ged$GWNO == c_id,]
	c_tmp <-  cshp2000[cshp2000$GWCODE == c_id,]
  conflict <- to_ppp(g_tmp,c_tmp)	
  conflict_data_list[[k]] <- conflict
  model <- to_ppm(conflict,c_tmp,1)
  conflict_model_list_1[[k]] <- model
  model <- to_ppm(conflict,c_tmp,2)
  conflict_model_list_2[[k]] <- model
  model <- to_ppm(conflict,c_tmp,3)
  conflict_model_list_3[[k]] <- model
  model <- to_ppm(conflict,c_tmp,4)
  conflict_model_list_4[[k]] <- model
  model <- to_ppm(conflict,c_tmp,5)
  conflict_model_list_5[[k]] <- model
  model <- to_ppm(conflict,c_tmp,6)
  conflict_model_list_6[[k]] <- model
  model <- to_ppm_star(conflict,c_tmp)
	conflict_model_list_star[[k]] <- model
  conflict_name_list[[k]]  <- as.character(c_tmp$ISONAME)
  cat(k,"/",length(relevant_countries$GWCODE),"\n")
	k <- k + 1
}
#Warnings indicate NAs in the covariate grids

#Generate in-sample predictions for optimized model
for (i in 1:length(conflict_data_list)){
  simulate_events(conflict_data_list[[i]],conflict_model_list_star[[i]],paste(conflict_name_list[[i]],"_model",sep=""),n_simulations,FALSE)
}

#####################################
#Generate a coefficient table for si#
#####################################

parameters <- matrix(nrow=length(conflict_model_list_6),ncol=length(conflict_model_list_6[[1]]$coef))
for (i in 1:length(conflict_model_list_6[[1]]$coef)){
  for(j in 1:length(conflict_model_list_6)){
    parameters[j,i] <- conflict_model_list_6[[j]]$coef[i]
  }
}

tmp <- data.frame(parameters)
tmp <- cbind(unlist(conflict_name_list), tmp)
names(tmp) <- c("country",names(conflict_model_list_6[[j]]$coef))
xtable(tmp)

##############################
#In-sample prediction results#
##############################
#Please note that running the code below will take some time! 

is_progress <- c()
results_is <- is_predict_plot(conflict_data_list,conflict_model_list_1,conflict_name_list,FALSE,n_simulations) 
is_progress <- c(is_progress,1)
results_is <- cbind(results_is,is_predict_plot(conflict_data_list,conflict_model_list_2,conflict_name_list,FALSE,n_simulations)$cv)
is_progress <- c(is_progress,2)
results_is <- cbind(results_is,is_predict_plot(conflict_data_list,conflict_model_list_3,conflict_name_list,FALSE,n_simulations)$cv)
is_progress <- c(is_progress,3)
results_is <- cbind(results_is,is_predict_plot(conflict_data_list,conflict_model_list_4,conflict_name_list,FALSE,n_simulations)$cv)
is_progress <- c(is_progress,4)
results_is <- cbind(results_is,is_predict_plot(conflict_data_list,conflict_model_list_5,conflict_name_list,FALSE,n_simulations)$cv)
is_progress <- c(is_progress,5)
results_is <- cbind(results_is,is_predict_plot(conflict_data_list,conflict_model_list_6,conflict_name_list,FALSE,n_simulations)$cv)
is_progress <- c(is_progress,6)
results_is <- cbind(results_is,is_predict_plot(conflict_data_list,conflict_model_list_star,conflict_name_list,FALSE,n_simulations)$cv)
is_progress <- c(is_progress,7)
results_is <- cbind(results_is,is_predict_plot(conflict_data_list,conflict_model_list_6,conflict_name_list,TRUE,n_simulations*2)$cv)
is_progress <- c(is_progress,8)

names(results_is) <- c("country","m1","m2","m3","m4","m5","m6","m7","random")

results_is[,2:8] <- round(results_is[,2:8],2)

results_is <- rbind(results_is,c(0,0,0,0,0,0,0,0))
results_is[length(results_is[,1]),1] <- "sum"
results_is[length(results_is[,1]),2] <- sum(results_is$m1)
results_is[length(results_is[,1]),3] <- sum(results_is$m2)
results_is[length(results_is[,1]),4] <- sum(results_is$m3)
results_is[length(results_is[,1]),5] <- sum(results_is$m4)
results_is[length(results_is[,1]),6] <- sum(results_is$m5)
results_is[length(results_is[,1]),7] <- sum(results_is$m6)
results_is[length(results_is[,1]),8] <- sum(results_is$m7)
results_is[length(results_is[,1]),9] <- sum(results_is$random)

results_is <- rbind(results_is,c(0,0,0,0,0,0,0,0))
results_is[length(results_is[,1]),1] <- "normalized"
results_is[length(results_is[,1]),2] <- sum(results_is$m1) / sum(results_is$random)
results_is[length(results_is[,1]),3] <- sum(results_is$m2) / sum(results_is$random)
results_is[length(results_is[,1]),4] <- sum(results_is$m3) / sum(results_is$random)
results_is[length(results_is[,1]),5] <- sum(results_is$m4) / sum(results_is$random)
results_is[length(results_is[,1]),6] <- sum(results_is$m5) / sum(results_is$random)
results_is[length(results_is[,1]),7] <- sum(results_is$m6) / sum(results_is$random)
results_is[length(results_is[,1]),8] <- sum(results_is$m7) / sum(results_is$random)
results_is[length(results_is[,1]),9] <- sum(results_is$random) / sum(results_is$random)

for(i in 2:length(results_is[1,])){
  results_is[,i] <- round(as.numeric(results_is[,i]),2)
}

xtable(results_is)

#####################################
#Cross-validation prediction results#
#####################################

progress <-c()
results_cv_var_sel <- k_fold_cross_plot(conflict_data_list,conflict_model_list_1,conflict_name_list,FALSE,n_simulations)
progress <- c(progress,1)
results_cv_var_sel <- cbind(results_cv_var_sel,k_fold_cross_plot(conflict_data_list,conflict_model_list_2,conflict_name_list,FALSE,n_simulations)$cv)
progress <- c(progress,2)
results_cv_var_sel <- cbind(results_cv_var_sel,k_fold_cross_plot(conflict_data_list,conflict_model_list_3,conflict_name_list,FALSE,n_simulations)$cv)
progress <- c(progress,3)
results_cv_var_sel <- cbind(results_cv_var_sel,k_fold_cross_plot(conflict_data_list,conflict_model_list_4,conflict_name_list,FALSE,n_simulations)$cv)
progress <- c(progress,4)
results_cv_var_sel <- cbind(results_cv_var_sel,k_fold_cross_plot(conflict_data_list,conflict_model_list_5,conflict_name_list,FALSE,n_simulations)$cv)
progress <- c(progress,5)
results_cv_var_sel <- cbind(results_cv_var_sel,k_fold_cross_plot(conflict_data_list,conflict_model_list_6,conflict_name_list,FALSE,n_simulations)$cv)
progress <- c(progress,6)
results_cv_var_sel <- cbind(results_cv_var_sel,k_fold_cross_plot(conflict_data_list,conflict_model_list_star,conflict_name_list,FALSE,n_simulations)$cv)
progress <- c(progress,7)
results_cv_var_sel <- cbind(results_cv_var_sel,k_fold_cross_plot(conflict_data_list,conflict_model_list_5,conflict_name_list,TRUE,n_simulations*2)$cv)
progress <- c(progress,8)

names(results_cv_var_sel) <- c("country","m1","m2","m3","m4","m5","m6","m7","random")
results_cv <- results_cv_var_sel
results_cv[,2:8] <- round(results_cv[,2:8],2)

results_cv <- rbind(results_cv_var_sel,c(0,0,0,0,0,0,0,0))
results_cv[length(results_cv[,1]),1] <- "sum"
results_cv[length(results_cv[,1]),2] <- sum(results_cv_var_sel$m1)
results_cv[length(results_cv[,1]),3] <- sum(results_cv_var_sel$m2)
results_cv[length(results_cv[,1]),4] <- sum(results_cv_var_sel$m3)
results_cv[length(results_cv[,1]),5] <- sum(results_cv_var_sel$m4)
results_cv[length(results_cv[,1]),6] <- sum(results_cv_var_sel$m5)
results_cv[length(results_cv[,1]),7] <- sum(results_cv_var_sel$m6)
results_cv[length(results_cv[,1]),8] <- sum(results_cv_var_sel$m7)
results_cv[length(results_cv[,1]),9] <- sum(results_cv_var_sel$random)

results_cv <- rbind(results_cv,c(0,0,0,0,0,0,0,0))
results_cv[length(results_cv[,1]),1] <- "normalized"
results_cv[length(results_cv[,1]),2] <- sum(results_cv_var_sel$m1) / sum(results_cv_var_sel$random)
results_cv[length(results_cv[,1]),3] <- sum(results_cv_var_sel$m2) / sum(results_cv_var_sel$random)
results_cv[length(results_cv[,1]),4] <- sum(results_cv_var_sel$m3) / sum(results_cv_var_sel$random)
results_cv[length(results_cv[,1]),5] <- sum(results_cv_var_sel$m4) / sum(results_cv_var_sel$random)
results_cv[length(results_cv[,1]),6] <- sum(results_cv_var_sel$m5) / sum(results_cv_var_sel$random)
results_cv[length(results_cv[,1]),7] <- sum(results_cv_var_sel$m6) / sum(results_cv_var_sel$random)
results_cv[length(results_cv[,1]),8] <- sum(results_cv_var_sel$m7) / sum(results_cv_var_sel$random)
results_cv[length(results_cv[,1]),9] <- sum(results_cv_var_sel$random) / sum(results_cv_var_sel$random)

for(i in 2:length(results_cv[1,])){
	results_cv[,i] <- round(as.numeric(results_cv[,i]),2)
}
xtable(results_cv)

####################
#Make case overview#
####################
gwnos <- unique(ged$GWNO)
data_cases <- as.data.frame(cbind(gwnos,rep(NA,length(gwnos)),rep(NA,length(gwnos)),rep(NA,length(gwnos))))
names(data_cases) <- c("gwno","country","min","max")
for (i in 1:length(data_cases[,1])){
	c_gwno <- data_cases$gwno[i]
	data_cases$country[i] <- unique(as.character(ged$COUNTRY[ged$GWNO ==  c_gwno]))
	data_cases$min[i] <- as.character(min(as.Date(ged$DATE_START)[ged$GWNO ==  c_gwno]))
	data_cases$max[i] <- as.character(max(as.Date(ged$DATE_START)[ged$GWNO ==  c_gwno]))
}
xtable(data_cases)

##################################
#Generate model averaging results#
##################################
#Please note that running the below code will take a long time! 

ma_res <- data.frame(matrix(nrow=6,ncol=12))
names(ma_res) <- c("variable","mean_cv",unlist(conflict_name_list))
ma_res$variable <- c("pop","lsg","acc","brd","nrd","veg")
for (k in 1:length(conflict_data_list)){
  #generating in-sample results
  tmp <- get_var_gain(conflict_model_list_6[[k]],conflict_data_list[[k]],TRUE,n_simulations)
  ma_res[1,k+2] <- mean(subset(tmp,pop == 1)$cv)
  ma_res[2,k+2] <- mean(subset(tmp,lsg == 1)$cv)
  ma_res[3,k+2] <- mean(subset(tmp,acc == 1)$cv)
  ma_res[4,k+2] <- mean(subset(tmp,brd == 1)$cv)
  ma_res[5,k+2] <- mean(subset(tmp,nrd == 1)$cv)
  ma_res[6,k+2] <- mean(subset(tmp,veg == 1)$cv)
}
for (k in 1:length(ma_res$variable)){
  ma_res$mean_cv[k] <- mean(as.numeric(ma_res[k,3:12]))
}
xtable(ma_res,digits=c(rep(4,13)))

